1 Introducción

Los datos que se van a analizar en este proyecto han sido obtenidos desde Kaggle. Contienen precios de casas que fueron vendidas desde mayo de 2014 hasta mayo de 2015 en King County que es un condado ubicado en el estado estadounidense de Washington.

2 Objetivo del estudio

Lo que queremos hacer con estos datos es clasificar/predecir las viviendas conforme a su precio dependiendo de las variables recogidas de cada una. Para ello se implementarán varios modelos, comprobando cuáles de ellos funciona mejor con los tipos de variables que contamos. Finalmente se evaluarán los modelos en base a distintas métricas.

3 Datos

3.1 Categorización del precio

En nuestro estudio inicial, la variable “Precio” es continua, por lo que vamos a categorizarla. Para decidir las categorizaciones se ha usado los cuantiles. Se van a realizar dos tipos de categorizaciones:

Categorización 1: se ha categorizado en dos grupos:
- B1: casas baratas (casas con un precio < 500.000).
- C1: casas caras (casas con un precio > 500.000).
Categorización 2: se ha categorizado en tres grupos:
- B2: casas baratas (casas con un precio < 330.000).
- M2: casas con precio medio (330.000 < precio < 650.000).
- C2: casas caras (precio > 650.000).

En el siguiente mapa se puede visualizar cómo se distribuyen las casas “Caras” y “Baratas”. Se observa que las casas que están cercanas al agua y cerca de Seattle por la parte norte son más caras y hacia el sur son más baratas.

center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red"), 
                       datos$price_categ1 )

leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(datos$price_categ1))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat, zoom = 9) %>%

  addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ1,
            title = "Tipos de Casas",
            opacity = 1)

En este otro mapa se puede visualizar cómo se distribuyen las casas según el precio en tres categorías:“Caras”, “Medio” y “Baratas”. Se puede observar cómo con esta categorización no está tan clara la separación entre clases.

center_lon = median(datos$long,na.rm = TRUE)
center_lat = median(datos$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red","yellow"), 
                       datos$price_categ2 )

leaflet(datos) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
  addCircles(lng = ~long, lat = ~lat,
             color = ~factpal2(datos$price_categ2))  %>%
  # controls
  setView(lng=center_lon, lat=center_lat,zoom = 9) %>%

  addLegend("bottomright", pal = factpal2 , values = ~datos$price_categ2,
            title = "Tipos de Casas 3 categorías",
            opacity = 1)

Por lo tanto, viendo ambos mapas, nos quedamos con la primera categorización: casas baratas y caras.

# Se elimina la segunda categorización
datos$price_categ2= NULL

3.2 Train, test y validación

Se va a separar los datos en los 3 conjuntos de datos fundamentales:

Conjunto de datos de entrenamiento: en nuestro estudio datos_train, se corresponde con el 70% del total de los datos.
Conjunto de datos de test: en nuestro estudio datos_test, se corresponde con el 15% del total de los datos.
Conjunto de datos de validación: en nuestro estudio datos_validacion, se corresponde con el 15% del total de los datos.

num_total=nrow(datos)
set.seed(122556) #reproductividad

# 70% para train
indices_train = sample(1:num_total, .7*num_total)
datos_train = datos[indices_train,]

# 15% para test
indices=seq(1:num_total)
indices_test=indices[-indices_train]
indices_test1 = sample(indices_test, .15*num_total)
datos_test = datos[indices_test1,]

# 15% para validacion
indices_validacion=indices[c(-indices_train,-indices_test1)]
datos_validacion=datos[indices_validacion,]

3.3 Análisis exploratorio

Se van a realizar transformaciones de un conjunto de variables, estas transformaciones se aplicarán a cada conjunto de datos: train, test y validación.

Se realiza una transformación logarítmica sobre las variables sqft_living (pies cuadrados de la casa), sqft_lot (pies cuadrados del jardín) y sqft_above (pies cuadrados por encima del suelo). Hay que aclarar que esta última variable es la diferencia entre sqft_living y sqft_basement, por lo que va a estar altamentente correlada con sqft_living.
Se categorizan las variables:
- Bathroom, esta varible puede tomar valores decimales de 0.25 en 0.25. El número de baños se contabiliza por las piezas y cada baño completo tiene 4 piezas. Por lo que con la nueva agrupación toma valores de 1 a 8 baños.
- Sqft_basement, se categoriza como 0 las casas que no tienen sótano y 1 las casas que sí tienen sótano.
- Grade, se va a categorizar del siguiente modo: con valor 0: calidad Baja, 1: calidad media y 2: calidad alta.
- Year_renovated, se categoriza como 0: no ha tenido renovación y 1: sí ha tenido renovación.
Se pasan a factor las variables: waterfront, view, condition, grade_categ y zipcode.
Se eliminan Outliers.

3.3.1 Transformaciones datos Train

Se realizan las transformaciones anteriormente mencionada y se eliminan outliers.

datos_train <- datos_train[,-2]

datos_train$id <- as.factor(datos_train$id)

datos_train$bathrooms_group <- cut(datos_train$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_train$bathrooms_group <- as.numeric(as.character(datos_train$bathrooms_group))

datos_train$log_sqft_living <- log10(datos_train$sqft_living)
datos_train$log_lot <- log10(datos_train$sqft_lot)
datos_train$log_above <- log10(datos_train$sqft_above)

datos_train$sqft_basement_cat <- cut(datos_train$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_train$waterfront<-as.factor(datos_train$waterfront)

datos_train$view<-as.factor(datos_train$view)

datos_train$condition<-as.factor(datos_train$condition)

datos_train$grade_categ <- cut(datos_train$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_train$yr_renovated_catg <-cut(datos_train$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_train$zipcode<-as.factor(datos_train$zipcode)

# Eliminación de Outliers

datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_habitaciones<-datos_train[datos_train$bedrooms==0,]$posicion
datos_train<-datos_train[-indices_cero_habitaciones,]

datos_train$posicion<-c(1:nrow(datos_train))
indices_cero_banos<-datos_train$posicion[datos_train$bathrooms_group==0]
datos_train<-datos_train[-indices_cero_banos,]

datos_train$posicion<-c(1:nrow(datos_train))
indice_hab33 <- datos_train[datos_train$bedrooms==33,]$posicion
datos_train[datos_train$posicion == indice_hab33,]$bedrooms = 3

Transformaciones adicionales:

Se han realizado dos categorizaciones adicionales sobre el conjunto de datos. Después de implementar varios modelos, se llegó a la conclusión de que algunas variables podían mejorar los resultados de los modelos siendo agrupadas. Para analizar cómo recategorizar estas variables se ha usado un árbol de decisión. Las variables son: zipcode y bathrooms_group.

Zipcode, esta variable es de tipo factor y tenía 70 códigos postales, por lo que se ha decidido aplicar un árbol de decisión para ver cómo clasificaba los códigos postales y así volver a categorizarla según el resultado obtenido.

set.seed(1234)
model_selec_zipcode<-rpart(price_categ1~zipcode,data=datos_train ,parms=list(split="gini"))
print(model_selec_zipcode)

## n= 15119 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 15119 6385 B1 (0.5776837 0.4223163)  
##   2) zipcode=98001,98002,98003,98010,98011,98014,98019,98022,98023,98024,98028,98030,98031,98032,98034,98038,98042,98045,98055,98056,98058,98059,98070,98092,98106,98108,98118,98125,98126,98133,98144,98146,98148,98155,98166,98168,98178,98188,98198 8243 1336 B1 (0.8379231 0.1620769) *
##   3) zipcode=98004,98005,98006,98007,98008,98027,98029,98033,98039,98040,98052,98053,98065,98072,98074,98075,98077,98102,98103,98105,98107,98109,98112,98115,98116,98117,98119,98122,98136,98177,98199 6876 1827 C1 (0.2657068 0.7342932) *

Se va a categorizar en dos: Zona1 y Zona2.

datos_train$zona<-recode(datos_train$zipcode, "98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_train$zipcode = NULL

bathrooms_group, se aplica el mismo método que con zipcode para ver cómo se puede categorizar esta variable. Toma valores de 1 a 8 y queremos reducir el número de niveles.

set.seed(1234)
model_selec_bathrooms<-rpart(price_categ1~bathrooms_group,data=datos_train ,parms=list(split="gini"))
# print(model_selec_bathrooms)
fancyRpartPlot(model_selec_bathrooms, caption='')

En el resultado del modelo se ve que corta en el número de baños < 2.5, por lo que se va a categorizar como 0 aquellas casas que tengan de 1 a 2 baños y como 1 las casas que tengan más de 2 baños.

datos_train$bathrooms_group <- cut(datos_train$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))

# Limpiamos el dataframe:
datos_train_limpio <- datos_train[c(3,21:25,8:10,26,16,17,29,27,20)]

#Eliminamos sqft_above:
datos_train_limpio$log_above = NULL

datos_train_numeric <- datos_train_limpio %>% select_if(is.numeric)

3.3.2 Transformaciones datos Test

Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de test.

datos_test <- datos_test[,-2]

datos_test$id <- as.factor(datos_test$id)

datos_test$bathrooms_group <- cut(datos_test$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_test$bathrooms_group <- as.numeric(as.character(datos_test$bathrooms_group))

datos_test$log_sqft_living <- log10(datos_test$sqft_living)
datos_test$log_lot <- log10(datos_test$sqft_lot)
datos_test$log_above <- log10(datos_test$sqft_above)

datos_test$sqft_basement_cat <- cut(datos_test$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_test$waterfront<-as.factor(datos_test$waterfront)

datos_test$view<-as.factor(datos_test$view)

datos_test$condition<-as.factor(datos_test$condition)

datos_test$grade_categ <- cut(datos_test$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_test$yr_renovated_catg <-cut(datos_test$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_test$zipcode<-as.factor(datos_test$zipcode)

#codificar la variable Zipcode

datos_test$zona<-recode(datos_test$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_test$zipcode = NULL
datos_test$bathrooms_group <- cut(datos_test$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
datos_test_limpio <- datos_test[c(3,21:25,8:10,26,16,17,28,27,20)]
datos_test_limpio$log_above = NULL
datos_test_numeric <- datos_test_limpio %>% select_if(is.numeric)

3.3.3 Transformaciones datos Validación

Realizamos todas las transformaciones y categorizaciones para el conjunto de datos de Validación.

datos_validacion <- datos_validacion[,-2]

datos_validacion$id <- as.factor(datos_validacion$id)

datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms,breaks = c(-1,0.25,1,2,3,4,5,6,7,8),labels=c(0,1,2,3,4,5,6,7,8))
datos_validacion$bathrooms_group <- as.numeric(as.character(datos_validacion$bathrooms_group))

datos_validacion$log_sqft_living <- log10(datos_validacion$sqft_living)
datos_validacion$log_lot <- log10(datos_validacion$sqft_lot)
datos_validacion$log_above <- log10(datos_validacion$sqft_above)

datos_validacion$sqft_basement_cat <- cut(datos_validacion$sqft_basement,breaks = c(-1,0,6000),labels=c(0,1))

datos_validacion$waterfront<-as.factor(datos_validacion$waterfront)

datos_validacion$view<-as.factor(datos_validacion$view)

datos_validacion$condition<-as.factor(datos_validacion$condition)

datos_validacion$grade_categ <- cut(datos_validacion$grade, breaks = c(0,4,9,13), labels = c(0,1,2))

datos_validacion$yr_renovated_catg <-cut(datos_validacion$yr_renovated, breaks=c(-0.5,1933, 2015), labels= c("0","1"))

datos_validacion$zipcode<-as.factor(datos_validacion$zipcode)

#codificar la variable Zipcode

datos_validacion$zona<-recode(datos_validacion$zipcode, " 98001=1; 98002=1; 98003=1; 98010=1; 98011=1; 98014=1; 98019=1; 98022=1; 98023=1; 98024=1; 98028=1; 98030=1; 98031=1; 98032=1; 98034=1; 98038=1; 98042=1; 98045=1; 98055=1; 98056=1; 98058=1; 98059=1; 98070=1; 98092=1 ;98106=1; 98108=1; 98118=1; 98125=1; 98126=1; 98133=1; 98144=1; 98146=1; 98148=1; 98155=1; 98166=1; 98168=1; 98178=1; 98188=1; 98198=1; 98004=2; 98005=2; 98006=2; 98007=2; 98008=2; 98027=2; 98029=2; 98033=2; 98039=2; 98040=2; 98052=2; 98053=2; 98065=2; 98072=2; 98074=2; 98075=2; 98077=2; 98102=2; 98103=2; 98105=2; 98107=2; 98109=2; 98112=2; 98115=2; 98116=2; 98117=2; 98119=2; 98122=2; 98136=2; 98177=2; 98199=2")

datos_validacion$zipcode = NULL
datos_validacion$bathrooms_group <- cut(datos_validacion$bathrooms_group, breaks = c(-1,2.5,8),labels=c(0,1))
datos_validacion_limpio <- datos_validacion[c(3,21:25,8:10,26,16,17,28,27,20)]
datos_validacion_limpio$log_above = NULL
datos_validacion_numeric <- datos_validacion_limpio %>% select_if(is.numeric)

3.3.4 Resumen EDA

La distribución de las casas según la nueva categorización daría como resultado un 57,8% de casas baratas (B1) y un 42,23% de casas caras (C1):

datos_train_limpio %>% 
  group_by(price_categ1) %>% 
  summarise(Count = n())%>% 
  mutate(percent = prop.table(Count)*100)%>%
  ggplot(aes(reorder(price_categ1, - percent), percent), fill = price_categ1) +
  geom_col(fill = c("salmon", "cyan3"), size = 0.5) +
  xlab("Precio de las casas") + 
  ylab("Porcentaje") +
  ggtitle("Porcentaje de casas por precio") +
  theme(plot.title = element_text(hjust = 0.5))

En cuanto a la relación de las variables que hacen referencia a las características de las casas con la variable respuesta categorizada, se observa lo siguiente:

p1<-ggplot(data = datos_train_limpio)+
  geom_bar(aes(x=bedrooms,fill=price_categ1,bins=30, position = "identity")) 
p2<-ggplot(data=datos_train_limpio)+
  geom_bar(aes(x=bathrooms_group,fill = price_categ1))
p3<-ggplot(data=datos_train_limpio)+
  geom_density(aes(x=log_sqft_living, fill=price_categ1))
p4<-ggplot(data=datos_train_limpio)+
  geom_density(aes(x=log_lot, fill=price_categ1))+
  facet_grid(price_categ1~., scales = 'free')

grid.arrange(p1, p2, p3, p4, nrow = 2)

Por último, para la relación entre las variables que expresan la calidad de las casas y el precio, se obtienen los siguientes resultados:

p5<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=waterfront,fill=price_categ1, stat="count"))
p6<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=view,fill=price_categ1, stat="count"))
p7<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=condition,fill=price_categ1, stat="count"))
p8<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=grade_categ,fill=price_categ1, stat="count"))
p9<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=yr_renovated_catg,fill=price_categ1, stat="count"))
p10<-ggplot(data = datos_train_limpio) + geom_bar(aes(x=zona,fill=price_categ1,bins=30, position = "identity"))

grid.arrange(p5,p6,p7,p8,p9,p10, nrow = 3)

4 Análisis de componentes principales: PCA

Es un método que permite simplificar la complejidad de espacios muestrales con muchas dimensiones conservando la mayor cantidad de información posible.

datospca<-datos_train_numeric
pca<-prcomp(datospca,scale=T)
plot(pca, main = 'Análisis de componentes principales', xlab= 'Componente')

summary(pca)

## Importance of components:
##                          PC1    PC2    PC3    PC4     PC5
## Standard deviation     1.414 1.0885 0.9290 0.7850 0.57947
## Proportion of Variance 0.400 0.2369 0.1726 0.1232 0.06716
## Cumulative Proportion  0.400 0.6370 0.8096 0.9328 1.00000

pca$rotation

##                        PC1        PC2         PC3          PC4         PC5
## bedrooms         0.5170755 -0.4142446  0.36490069  0.108832991  0.64500949
## log_sqft_living  0.5825593 -0.3343737  0.09165417  0.008721066 -0.73507980
## log_lot          0.4612158  0.3615021 -0.29757928 -0.737528449  0.15522444
## lat             -0.1072636 -0.6381109 -0.75248337 -0.052915963  0.11080474
## long             0.4111352  0.4227604 -0.45128965  0.664327488  0.08513595

biplot(x = pca, scale = 0, cex = 0.6, col = c("blue4", "brown3"))

Se observa, como la componente 1 está dando un peso de 0.51 a bedrooms, 0.58 a log_sqft_living, 0.46 a log_lot, -0.10 a lat y 0.41 a long. La componente 2 -0.41 a bedrooms, -0.33 a log_sqft_living, 0.36 a log_lot, -0.63 a lat y 0.42 a long. Esta componente estaría diferenciando bedrooms, log_sqft_living y lat, frente a log_lot y long. Estas conclusiones se extraen de los pesos que cada autovector da a cada variable. Con las dos primeras componenetes se estarían explicando el 64% de la variabilidad de los datos. Si se consideraran las tres primeras dimensiones, se llegaría al 81%.

5 Métodos de agrupamiento-No supervisado

A continuación, se implementan tres métodos de agrupamiento no supervisado. Primero, un método jerárquico para averiguar (si es posible) la cantidad adecuada de clústers y, a continuación, K-Means y K-Medoids. La agrupación es una técnica para agrupar puntos de datos similares en un grupo y separar las diferentes observaciones en diferentes grupos, de forma que un mismo cluster agrupe casas con características homogéneas que se diferencien en cierta medida del resto de clústers.

5.1 Clúster Jerárquico

En el Clustering Jerárquico los clústers se crean de manera que tengan un orden predeterminado. En nuestro estudio se va a aplicar un método aglomerativo que consiste en que cada observación se asigna a su propio clúster. Luego, se calcula la similitud (o distancia) entre cada uno de los clústers y los dos clústers más similares se fusionan en uno. Finalmente, los pasos 2 y 3 se repiten hasta que solo quede un grupo.

Aplicando este método a nuestros datos queremos observar si siguen algún patrón para poder agrupar las casas.

Escalamos las variables numéricas, es decir, cada variable ahora tendrá una media cero y una desviación estándar uno. Por otro lado, se halla la matriz de distancias mediante la función dist que usa por defecto la ‘Euclidea’.

datos_scale<- as.data.frame(scale(datos_train_numeric))
matriz_dist=dist(datos_scale)

En este caso particular, probamos con dos clústers dado que la categorización del precio se hace en base a dos clases (caras y baratas):

set.seed(1234)
modelo_hc1= hclust(matriz_dist, method = "average")
plot(modelo_hc1, sub='')
rect.hclust(modelo_hc1, k=2, border="red")

set.seed(1234)
modelo_hc2= hclust(matriz_dist, method = "average")
plot(modelo_hc2, sub='')
rect.hclust(modelo_hc2, k=10, border="red")

A continuación vemos cómo se han agrupado los datos marcando que el número de clústers sea 2. Como se puede ver, prácticamente todas las casas están en el grupo 1. Si se considera $k=10$, sucede lo mismo, se agrupan prácticamente el total de observaciones en el grupo 1 y prácticamente ninguna observación en el resto de grupos.

grupos2=cutree(modelo_hc1,k=2)
table(datos_train_limpio$price_categ1, grupos2)

##     grupos2
##         1    2
##   B1 8724   10
##   C1 6385    0

grupos10=cutree(modelo_hc2,k=10)
table(datos_train_limpio$price_categ1, grupos10)

##     grupos10
##         1    2    3    4    5    6    7    8    9   10
##   B1 8646    8   43    2   10    4    0   18    0    3
##   C1 6288   32   37    7    0   19    1    0    1    0

Se concluye que la agrupación en dos clases está muy descompensada, obteniendo un total de 15.109 casas en el grupo 1 y 10 en el grupo 2 y lo mismo sucede con 10 grupos. Dados estos resultados, se decide que no es un método adecuado para los datos con los que contamos y por tanto, no se evaluan los conjuntos de test y validación.

5.2 Clúster no Jerárquico- K-Means

El método de K-Means basa su funcionamiento en agrupar los datos de entrada en un total de k conjuntos definidos por un centroide, cuya distancia con los puntos que pertenecen a cada uno de los datos sea la menor posible.

Primero se van a realizar los gráficos para ver cómo están diferenciadas las casas por precio. Para poder visualizarlo en dos dimensiones se ha usado la función “prcomp” (PCA).

colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]

plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n",main= "Dos categorías")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(datos_train_limpio$price_categ1), col=colores11)

A continuación se va a aplicar el método de agrupamiento K-means, al igual que en el método anterior se le pasa la matriz de distacias y se van a agrupar los datos en dos conjuntos:

set.seed(1234)
model_km <- kmeans(matriz_dist, centers=2)

Se representa mediante PCA el resultado obtenido con $k=2$, y se observa como en ambos grupos ha incluido tanto casas caras como baratas.

colores1= c("red","blue")
colores11 = colores1[datos_train_limpio$price_categ1]

plot(prcomp(datos_train_numeric, scale = T)$x[,1:2], type="n")
text(prcomp(datos_train_numeric, scale = T)$x[,1:2], as.character(model_km$cluster), col=colores11)

Cálculo del k-óptimo.

Se va a determinar la cantidad óptima de centroides a utilizar a partir del Método del Codo. Para ello, aplicaremos la función kmeans al conjunto de datos, variando en cada caso el valor de k y acumulando los valores de WCSS (Within-Cluster-Sum-of-Squares) obtenidos:

set.seed(1234)
wcss <- vector()
for(i in 1:20){
wcss[i] <- sum(kmeans(datos_scale, i)$withinss)
}

ggplot() + geom_point(aes(x = 1:20, y = wcss), color = 'blue') + 
geom_line(aes(x = 1:20, y = wcss), color = 'blue') + 
ggtitle("Método del Codo") + 
xlab('Cantidad de Centroides k') + 
ylab('WCSS')

Como se observa en la gráfica, el K-óptimo que se podría aplicar sería de 10. Después de volver a implementar el modelo con $k=10$, se concluye que sigue sin haber un patrón que permita agrupar las casas en función de sus características.

Se comparan los datos obtenido con $k=2$ y con $k=10$:

load("model_km10.RData")
table(model_km$cluster, datos_train_limpio$price_categ1) #asignación de observación a los cluster

##    
##       B1   C1
##   1 1554 2372
##   2 7180 4013

table(model_km10$cluster, datos_train_limpio$price_categ1)

##     
##        B1   C1
##   1   141  249
##   2   146  862
##   3   236 1500
##   4  1209  107
##   5  1545 1053
##   6  1243  808
##   7  1587   35
##   8   610  189
##   9  1493  315
##   10  524 1267

A continuación, para visualizar si el agrupamiento que se ha llevado a cabo está relacionado con el precio de las casas (Caras, Baratas), se representa en el mapa (para el caso de $k=2$) que se muestra a continuación.

clusterkmeans=as.data.frame(model_km$cluster)
clusterkmeans$indice= as.integer(rownames(clusterkmeans))

colnames(clusterkmeans)[1]= "categoria_price_km"
clustering1= clusterkmeans[order(clusterkmeans$indice),]


center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)

factpal1 <- colorFactor(c("green","red"),
                    clustering1$categoria_price_km )

leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
         color = ~factpal1(clustering1$categoria_price_km))  %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 9) %>%

addLegend("bottomright", pal = factpal1 , values = ~clustering1$categoria_price_km,
        title = "Tipos de casas",
        opacity = 1)

Comparando el mapa con el de los datos iniciales, dista bastante. Por lo que deducimos que la agrupación que está realizando K-Means con $k=2$ y con $k=10$ no sigue ningún patrón. Como ya se suponía, la naturaleza de los datos no permiten una agrupación de los mismos, por ejemplo, una casa cara puede tener el mismo número de pies cuadrado que una barata, que esté situada en otro barrio.

5.3 K-Medoids

K-medoids es un método de clustering muy similar a K-means en cuanto a que ambos agrupan las observaciones en K clústers, donde K es un valor preestablecido. La diferencia es que, en K-medoids, cada clúster está representado por una observación presente en el clúster (medoid). En nuestro estudio será una observación de una casa, mientras que en K-means cada clúster está representado por su centroide, que se corresponde con el promedio de todas las observaciones del clúster pero con ninguna casa en particular.

Este algoritmo es menos sensible al ruido y a los valores atípicos en comparación con k-means, porque usa medoides como centros de clúster en lugar de centroides (utilizados en k-means).

En primer lugar probamos con las variables (bedrooms, log_sqft_living, log_lot, lat, long y zona), y en segundo lugar incluimos todas las variables para comparar:

datoskmedoids1 = datos_train_limpio[,c(1,3,4,10:12)]
model_medoids1 = pam(x = datoskmedoids1, k = 2, keep.diss = TRUE, keep.data = TRUE)
model_medoids1$medoids

##      bedrooms log_sqft_living  log_lot     lat     long zona
## 3470        3        3.209515 3.854306 47.5048 -122.241    1
## 8813        4        3.399674 3.998390 47.5973 -122.177    2

datoskmedoids2 = datos_train_limpio[,-14]
model_medoids2 = pam(x = datoskmedoids2, k = 2, keep.diss = TRUE, keep.data = TRUE)
model_medoids2$medoids

##       bedrooms bathrooms_group log_sqft_living  log_lot sqft_basement_cat
## 9606         3               1        3.152288 3.870404                 1
## 18626        4               2        3.411620 3.865933                 1
##       waterfront view condition grade_categ     lat     long zona
## 9606           1    1         3           2 47.4949 -122.239    1
## 18626          1    1         3           2 47.5647 -122.090    2
##       yr_renovated_catg
## 9606                  1
## 18626                 1

Una vez más se representa el mapa para comprobar los resultados de la agrupación y pese a que los medoids aparecían separados, la clasificación es muy heterogénea y no se corresponde con la categorización del precio, por tanto se descarta este modelo.

clustering = sort(model_medoids1$clustering)
clustering = as.data.frame(model_medoids1$clustering)
clustering$indice= as.integer(rownames(clustering))

colnames(clustering)[1]= "categoria_price"
clustering2 = clustering[order(clustering$indice),]

center_lon = median(datoskmedoids1$long,na.rm = TRUE)
center_lat = median(datoskmedoids1$lat,na.rm = TRUE)

factpal2 <- colorFactor(c("green","red"), 
                      clustering2$categoria_price )

leaflet(datoskmedoids1) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
           color = ~factpal2(clustering2$categoria_price))  %>%
# controls
setView(lng=center_lon, lat=center_lat,zoom = 9) %>%

addLegend("bottomright", pal = factpal2 , values = ~clustering2$categoria_price,
          title = "Tipos de casas",
          opacity = 1)

6 Aprendizaje supervisado

6.1 GLM-Regresión Logística.

Con este modelo se va a estudiar si existe relación entre el hecho de que una casa sea “cara” ó “barata” dependiendo de las características de las casas. Se va a generar un modelo en el que a partir de las variables prediga la probabilidad de que una casa sea barata o cara.

6.1.1 Train

Inicialmente, se probó con todas las variables y después de implementar y evaluar varios modelos, se concluyó que las variables condition y grade_categ no eran importantes, por tanto se eliminaron para este modelo.

datos_train_rl <- datos_train_limpio[,-c(8,9)]
datos_train_rl$price_categ1<- recode(datos_train_rl$price_categ1, "'B1'=0; 'C1'=1")

Se seleccionan las variables log_sqft_living y zona, en función de la importancia que el Random Forest asigna a cada variable.

set.seed(1234)
model_glm = glm(price_categ1 ~log_sqft_living+zona, family = binomial, data =datos_train_rl)
summary(model_glm)

## 
## Call:
## glm(formula = price_categ1 ~ log_sqft_living + zona, family = binomial, 
##     data = datos_train_rl)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.2060  -0.4965  -0.1352   0.4116   3.3879  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -38.63064    0.72115  -53.57   <2e-16 ***
## log_sqft_living  11.11677    0.21336   52.10   <2e-16 ***
## zona2             3.43940    0.05977   57.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 20593  on 15118  degrees of freedom
## Residual deviance: 10215  on 15116  degrees of freedom
## AIC: 10221
## 
## Number of Fisher Scoring iterations: 6

La métrica que se ha considerado para evaluar los modelos es la F1-medida, por un lado para las casas caras, y por otro, para las baratas. Después se realiza la media:

\[\frac{\text{F1-medida }_{caras} + \text{F1-medida }_{baratas}}{2}\]

predicciones <- ifelse(test = model_glm$fitted.values > 0.5, yes = 1, no = 0)

tabla_glm_train <- table(model_glm$model$price_categ1, predicciones, dnn = c("observaciones", "predicciones"))
tabla_glm_train

##              predicciones
## observaciones    0    1
##             0 7729 1005
##             1 1279 5106

#caras
pred_caras_glm_train=tabla_glm_train[2,2]/(tabla_glm_train[2,2]+tabla_glm_train[1,2])
rec_caras_glm_train=tabla_glm_train[2,2]/(tabla_glm_train[2,2]+tabla_glm_train[2,1])
F_caras_reglog=(2*pred_caras_glm_train*rec_caras_glm_train)/(pred_caras_glm_train+rec_caras_glm_train)
cat(c('F1 casas caras: ', F_caras_reglog), '\n')

## F1 casas caras:  0.817221510883483

#baratas
pred_baratas_glm_train=tabla_glm_train[1,1]/(tabla_glm_train[1,1]+tabla_glm_train[2,1])
rec_baratas_glm_train=tabla_glm_train[1,1]/(tabla_glm_train[1,1]+tabla_glm_train[1,2])
F_baratas_reglog=(2*pred_baratas_glm_train*rec_baratas_glm_train)/(pred_baratas_glm_train+rec_baratas_glm_train)
cat(c('F1 casas baratas: ', F_baratas_reglog), '\n')

## F1 casas baratas:  0.871265922669372

#F-Medida
F_reglog_train= (F_caras_reglog+F_baratas_reglog)/2
cat(c('F1 global: ', F_reglog_train), '\n')

## F1 global:  0.844243716776427

accuracy_reglog_train = (tabla_glm_train[1,1]+tabla_glm_train[2,2]) / (tabla_glm_train[1,1]+tabla_glm_train[1,2]+tabla_glm_train[2,1]+tabla_glm_train[2,2])
cat(c('Accuracy: ', accuracy_reglog_train), '\n')

## Accuracy:  0.848931807659237

El resultado de la F1-medida para el modelo GLM considerando las variables: log_sqft_living y zona, es de 0.84, y lo mismo para la accuracy. Este modelo, al ser supervisado, mejora bastante la clasificación de las casas en función de sus características.

6.2 GAM

Los modelos GAM (Modelos Aditivos Generalizados) son una extensión de los modelos lineales que permiten obtener ajustes no lineales empleando múltiples predictores. Los modelos lineales presentan la restricción de que la predicción es una función lineal de los parámetros del modelo. En cambio, los modelos GAM usan funciones suavizadoras para ajustar los datos a funciones no paramétricas.

6.2.1 Train

set.seed(1234)

model_gam_categ <- gam(price_categ1 ~ 
                     s(log_sqft_living, bs = "ps") + 
                     s(lat, bs = "ps") + 
                     s(long, bs = "ps"), data = datos_train_limpio, binomial(link = "logit"))



summary(model_gam_categ)

## 
## Family: binomial 
## Link function: logit 
## 
## Formula:
## price_categ1 ~ s(log_sqft_living, bs = "ps") + s(lat, bs = "ps") + 
##     s(long, bs = "ps")
## 
## Parametric coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.73608    0.06029  -12.21   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                      edf Ref.df Chi.sq p-value    
## s(log_sqft_living) 5.832  5.988 2702.0  <2e-16 ***
## s(lat)             8.882  8.986 2827.9  <2e-16 ***
## s(long)            6.660  7.211  203.3  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.624   Deviance explained =   57%
## UBRE = -0.41153  Scale est. = 1         n = 15119

visreg(model_gam_categ)

Se evalúa este modelo mediante el cálculo de la diferencia entre el precio predicho y el observado. Para el caso de train, la diferencia media entre el precio real y el predicho es de 118.300$; para test es de 120.500$ y para valicación 116.000$.

predicciones <- predict(model_gam_categ, datos_train_limpio)
predicciones <- ifelse(test = predicciones > 0.5, yes = 1, no = 0)

tabla_gam_train <- table(model_gam_categ$model$price_categ1, predicciones, dnn = c("observaciones", "predicciones"))
tabla_gam_train

##              predicciones
## observaciones    0    1
##            B1 8140  594
##            C1 1526 4859

#caras
pred_caras_gam_train=tabla_gam_train[2,2]/(tabla_gam_train[2,2]+tabla_gam_train[1,2])
rec_caras_gam_train=tabla_gam_train[2,2]/(tabla_gam_train[2,2]+tabla_gam_train[2,1])
F_caras_gam=(2*pred_caras_gam_train*rec_caras_gam_train)/(pred_caras_gam_train+rec_caras_gam_train)
cat(c('F1 casas caras: ', F_caras_gam), '\n')

## F1 casas caras:  0.820915695218787

#baratas
pred_baratas_gam_train=tabla_gam_train[1,1]/(tabla_gam_train[1,1]+tabla_gam_train[2,1])
rec_baratas_gam_train=tabla_gam_train[1,1]/(tabla_gam_train[1,1]+tabla_gam_train[1,2])
F_baratas_gam=(2*pred_baratas_gam_train*rec_baratas_gam_train)/(pred_baratas_gam_train+rec_baratas_gam_train)
cat(c('F1 casas baratas: ', F_baratas_gam), '\n')

## F1 casas baratas:  0.884782608695652

#F-Medida
F_gam_train= (F_caras_gam+F_baratas_gam)/2
cat(c('F1 global: ', F_gam_train), '\n')

## F1 global:  0.85284915195722

accuracy_gam_train = (tabla_gam_train[1,1]+tabla_gam_train[2,2]) / (tabla_gam_train[1,1]+tabla_gam_train[1,2]+tabla_gam_train[2,1]+tabla_gam_train[2,2])
cat(c('Accuracy: ', accuracy_gam_train), '\n')

## Accuracy:  0.859779085918381

6.3 KNN

6.3.1 Train

Ajuste de Hiperparámetros.

Este método se basa en clasificar cada dato en un grupo en función de sus k vecinos más cercanos. Para ello se usa la distancia de cada nuevo punto a los ya existentes. Para ajustar y comparar este modelo, se han usado tres métodos:

Usando train.kknn.
Usando tune.
Usando diferentes k de 1 a 50 y comprobando los resultados de la F1-medida.

En primer lugar con la función train.kknn se obtiene de manera automática el mejor valor de k hasta un máximo de $k=30$.

set.seed(1234)

knn_1 <- train.kknn(price_categ1 ~ bedrooms+log_sqft_living+lat+long, data = datos_train_limpio, kmax = 30)
knn_1

## 
## Call:
## train.kknn(formula = price_categ1 ~ bedrooms + log_sqft_living +     lat + long, data = datos_train_limpio, kmax = 30)
## 
## Type of response variable: nominal
## Minimal misclassification: 0.1107216
## Best kernel: optimal
## Best k: 20

En segundo lugar se usa la función tune que itera hasta $k=30$ y muestra el error y la dispersión para cada valor de k.

set.seed(1234)

knn_2 <- tune.knn(x = scale(datos_train_numeric[,c(1,2,4,5)]),
                y = as.factor(datos_train_limpio$price_categ1),
                k = 1:30,
                tunecontrol = tune.control(sampling = "boot"))

knn_2$best.parameters

##     k
## 23 23

plot(knn_2, main = "Mejor k según tune")

En tercer lugar, se comprueba manualmente los resultados de la F1-medida con diferentes modelos desde $k=1$ hasta el $k=50$ y se representa el resultado.

set.seed(1234)

k_maximo=50

rango=1:k_maximo
f1_modelos=c()

for (i in rango){
  knn_3=knn.cv(scale(datos_train_numeric[,c(1,2,4,5)]), cl=as.factor(datos_train_limpio$price_categ1),k=i)
  tabla=table(datos_train_limpio$price_categ1,knn_3)
  # f1 casas caras
  pred_means_caras=tabla[2,2]/(tabla[2,2]+tabla[1,2])
  rec_means_caras=tabla[2,2]/(tabla[2,2]+tabla[2,1])
  f1_caras=(2*pred_means_caras*rec_means_caras)/(pred_means_caras+rec_means_caras)
  # f1 casas baratas
  pred_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[2,1])
  rec_means_baratas=tabla[1,1]/(tabla[1,1]+tabla[1,2])
  f1_baratas=(2*pred_means_baratas*rec_means_baratas)/(pred_means_baratas+rec_means_baratas)
  f1_total = (f1_baratas + f1_caras)/2
  f1_modelos=c(f1_modelos,f1_total)
}


plot(f1_modelos)

cat(c('Valor óptimo de k: ', which.max(f1_modelos)))

## Valor óptimo de k:  14

De la primer forma obtenemos un valor óptimo de $k=20$, del segundo modo a partir de aproximadamente $k=23$ el error es muy pequeño, y de la forma manual, concluimos que el mejor valor de k en función de la F1-medida es de $k=14$. Finalmente, elegimos el valor de $k=14$, implementamos el modelo y lo evaluamos.

model_knn=knn.cv(scale(datos_train_numeric[,c(1,2,4,5)]), cl=as.factor(datos_train_limpio$price_categ1), k=14, prob = TRUE)

center_lon = median(datos_train_limpio$long,na.rm = TRUE)
center_lat = median(datos_train_limpio$lat,na.rm = TRUE)

factpal1 <- colorFactor(c("green","red"), 
                      model_knn)

leaflet(datos_train_limpio) %>% addProviderTiles("Esri.NatGeoWorldMap") %>%
addCircles(lng = ~long, lat = ~lat,
           color = ~factpal1(model_knn))  %>%
# controls
setView(lng=center_lon, lat=center_lat, zoom = 9) %>%

addLegend("bottomright", pal = factpal1 , values = ~model_knn,
          title = "Precio (en miles de $)",
          opacity = 1)

tabla_knn_train=table(datos_train_limpio$price_categ1, model_knn)
tabla_knn_train

##     model_knn
##        B1   C1
##   B1 7917  817
##   C1  914 5471

pred_caras_knn_train = tabla_knn_train[2,2]/(tabla_knn_train[2,2]+tabla_knn_train[1,2])
rec_caras_knn_train = tabla_knn_train[2,2]/(tabla_knn_train[2,2]+tabla_knn_train[2,1])
F_caras_knn_train=(2*pred_caras_knn_train*rec_caras_knn_train)/(pred_caras_knn_train+rec_caras_knn_train)
cat(c('F1 caras: ', F_caras_knn_train), '\n')

## F1 caras:  0.863410400063126

pred_baratas_knn_train=tabla_knn_train[1,1]/(tabla_knn_train[1,1]+tabla_knn_train[2,1])
rec_baratas_knn_train=tabla_knn_train[1,1]/(tabla_knn_train[1,1]+tabla_knn_train[1,2])
F_baratas_knn_train=(2*pred_baratas_knn_train*rec_baratas_knn_train)/(pred_baratas_knn_train+rec_baratas_knn_train)
cat(c('F1 baratas: ', F_baratas_knn_train), '\n')

## F1 baratas:  0.901451750640478

F_knn_train= (F_caras_knn_train+F_baratas_knn_train)/2
cat(c('F1 global: ', F_knn_train), '\n')

## F1 global:  0.882431075351802

accuracy_knn_train= (tabla_knn_train[1,1]+tabla_knn_train[2,2]) / (tabla_knn_train[1,1]+tabla_knn_train[1,2]+tabla_knn_train[2,1]+tabla_knn_train[2,2])
cat(c('Accuracy: ', accuracy_knn_train), '\n')

## Accuracy:  0.885508300813546

Como los valores de k óptimos son muy parecidos entre las diferentes técnicas usadas para determinar este parámetro, se elige $k=14$, con los que se obtiene una F1-medida y accuracy de 0.88.

6.4 Decision Tree

Los árboles de decisión son un método usado en distintas disciplinas como modelo de predicción. Este método es similar a diagramas de flujo, en los que llegamos a puntos donde se toman decisiones de acuerdo a una regla. De manera general, lo que hace este algoritmo es encontrar la variable independiente que mejor separa nuestros datos en grupos, que corresponden con las categorías de la variable objetivo, en nuestro caso casas caras frente a baratas.

6.4.1 Train

Como en este modelo se toman decisiones en base a puntos de corte en las diferentes variables, se ha decidido eliminar la latitud y longitud dado que ya tenemos una variable zona que aporta esta misma información.

El parámetro cp (parámetro de complejidad) da una idea de lo que se penaliza añadir una nueva rama en el modelo, un cp más pequeño implica que el tamaño del árbol sea más reducido y viceversa. Cualquier división que no mejore este parámetro se eliminará mediante validación cruzada. El parámetro minbucket indica el número de casas que tienen que quedar en cualquier nodo terminal del árbol. Por último, maxdeph indica la profundidad del árbol, o número máximo de capas, contando la raiz como profundidad = 0. Se categoriza para este modelo log_sqft_living que ahora toma los niveles P (<3.4) y G (>3.4) y se consideran las variables log_sqft_categ, zona y log_lot.

set.seed(1234)

datos_decision_tree<-datos_train_limpio
datos_decision_tree$log_sqft_categ <- cut(datos_decision_tree$log_sqft_living, breaks = c(0, 3.4, 5), labels = c("P", "G"))

decisiontree_model=rpart(price_categ1~log_sqft_categ+zona+log_lot, 
                       data=datos_decision_tree, 
                       parms=list(split="gini"),
                       control = rpart.control(cp = 0.005, maxdepth = 4, minbucket = 400))

# print(decisiontree_model)

fancyRpartPlot(decisiontree_model, caption='')

A continuación, se evalúa este modelo y el resultado para la F1-medida y accuracy es de 0.81.

tabla_train_arbol=table(obs = datos_decision_tree$price_categ1, pred = predict(decisiontree_model, type = "class"))
tabla_train_arbol

##     pred
## obs    B1   C1
##   B1 7045 1689
##   C1 1149 5236

pred_caras_dt_train = tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[1,2])
rec_caras_dt_train = tabla_train_arbol[2,2]/(tabla_train_arbol[2,2]+tabla_train_arbol[2,1])
F_caras_dt_train=(2*pred_caras_dt_train*rec_caras_dt_train)/(pred_caras_dt_train+rec_caras_dt_train)
cat(c('F1 caras: ', F_caras_dt_train), '\n')

## F1 caras:  0.786776859504132

pred_baratas_dt_train = tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[2,1])
rec_baratas_dt_train = tabla_train_arbol[1,1]/(tabla_train_arbol[1,1]+tabla_train_arbol[1,2])
F_baratas_dt_train=(2*pred_baratas_dt_train*rec_baratas_dt_train)/(pred_baratas_dt_train+rec_baratas_dt_train)
cat(c('F1 baratas: ', F_baratas_dt_train), '\n')

## F1 baratas:  0.832348771266541

F_dt_train= (F_caras_dt_train+F_baratas_dt_train)/2
cat(c('F1 global: ', F_dt_train), '\n')

## F1 global:  0.809562815385336

accuracy_dt_train = (tabla_train_arbol[1,1]+tabla_train_arbol[2,2]) / (tabla_train_arbol[1,1]+tabla_train_arbol[1,2]+tabla_train_arbol[2,1]+tabla_train_arbol[2,2])
cat(c('Accuracy: ', F_dt_train), '\n')

## Accuracy:  0.809562815385336

6.5 Random Forest

Los Random forest son una combinación de árboles predictores tal que cada árbol depende de los valores de un vector aleatorio probado independientemente y con la misma distribución para cada uno de estos. Es una modificación sustancial de bagging que construye una larga colección de árboles no correlacionados y luego los promedia.

6.5.1 Train

La implementación de este modelo se realiza mediante la función randomForest, y los parámetros de esta función son: ntree que indica el número de árboles que se van a considerar; importance que indica si se considerará la importancia de los predictores; proximity que indica si se calculará o no la proximidad entre las filas (vectores de características de las casas); mtry para indicar el número de variables que serán seleccionadas aleatoriamente en cada división; y replace para indicar si las muestras serán tomadas con o sin reemplazamiento.

set.seed(1234)
randomforest_model_all=randomForest(price_categ1~., 
                                  data=datos_train_limpio, 
                                  parms=list(split="gini"),
                                  ntree=50, 
                                  importance = FALSE, 
                                  proximity = FALSE, 
                                  mtry=6, 
                                  replace = TRUE)
randomforest_model_all$importance

##                   MeanDecreaseGini
## bedrooms                 196.59064
## bathrooms_group          177.27454
## log_sqft_living         2085.83268
## log_lot                  582.63306
## sqft_basement_cat         68.87435
## waterfront                16.37023
## view                     136.30240
## condition                113.35641
## grade_categ               87.35018
## lat                     1485.34365
## long                     704.32839
## zona                    1671.25307
## yr_renovated_catg         24.18168

# print(randomforest_model_all)

set.seed(1234)
randomforest_model=randomForest(price_categ1~log_sqft_living+log_lot+lat+long, 
                              data=datos_train_limpio, 
                              parms=list(split="gini"),
                              ntree=50, 
                              importance = FALSE, 
                              proximity = FALSE, 
                              mtry=3, 
                              replace = TRUE)
print(randomforest_model)

## 
## Call:
##  randomForest(formula = price_categ1 ~ log_sqft_living + log_lot +      lat + long, data = datos_train_limpio, parms = list(split = "gini"),      ntree = 50, importance = FALSE, proximity = FALSE, mtry = 3,      replace = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 50
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 10.54%
## Confusion matrix:
##      B1   C1 class.error
## B1 7930  804  0.09205404
## C1  790 5595  0.12372749

Después de implementar el modelo, se evalua, obteniendo una F1-medida y accuracy de 0.9, el mejor resultado obtenido hasta el momento.

tabla_train_randomforest = table(obs = datos_train_limpio$price_categ1, pred = predict(randomforest_model, type = "class") )
tabla_train_randomforest

##     pred
## obs    B1   C1
##   B1 7930  804
##   C1  790 5595

mosaic(tabla_train_randomforest, shade = T, colorize = T, main="Matriz Confusión",gp = gpar(fill = matrix(c("cyan3", "salmon", "salmon", "cyan3"), 2, 2)))

pred_caras_rf_train = tabla_train_randomforest[2,2]/(tabla_train_randomforest[2,2]+tabla_train_randomforest[1,2])
rec_caras_rf_train = tabla_train_randomforest[2,2]/(tabla_train_randomforest[2,2]+tabla_train_randomforest[2,1])
F_caras_rf_train=(2*pred_caras_rf_train*rec_caras_rf_train)/(pred_caras_rf_train+rec_caras_rf_train)
cat(c('F1 caras: ', F_caras_rf_train), '\n')

## F1 caras:  0.875312891113892

pred_baratas_rf_train = tabla_train_randomforest[1,1]/(tabla_train_randomforest[1,1]+tabla_train_randomforest[2,1])
rec_baratas_rf_train=tabla_train_randomforest[1,1]/(tabla_train_randomforest[1,1]+tabla_train_randomforest[1,2])
F_baratas_rf_train=(2*pred_baratas_rf_train*rec_baratas_rf_train)/(pred_baratas_rf_train+rec_baratas_rf_train)
cat(c('F1 baratas: ', F_baratas_rf_train), '\n')

## F1 baratas:  0.908674229403002

F_rf_train= (F_caras_rf_train+F_baratas_rf_train)/2
cat(c('F1 global: ', F_rf_train), '\n')

## F1 global:  0.891993560258447

accuracy_rf_train = (tabla_train_randomforest[1,1]+tabla_train_randomforest[2,2]) / (tabla_train_randomforest[1,1]+tabla_train_randomforest[1,2]+tabla_train_randomforest[2,1]+tabla_train_randomforest[2,2])
cat(c('Accuracy: ', accuracy_rf_train), '\n')

## Accuracy:  0.894569746676367

6.6 SVM

Una máquina de vectores de soporte (SVM) es un algoritmo de aprendizaje supervisado. Este algoritmo construye un hiperplano óptimo en forma de superficie de decisión, de modo que el margen de separación entre las dos clases en los datos sea lo más amplia posible. Los vectores de soporte hacen referencia a un pequeño subconjunto de las observaciones de entrenamiento que se utilizan como soporte para la ubicación óptima de la superficie de decisión.

La función svm de la librería e1071, toma como parámatros el kernel a emplear (en nuestro caso se usa un kernel lineal, linear), el coste (cost) para ajustar el número de vectores soporte y probability para indicar si el modelo debe permitir predicciones de probabilidad.

Ajuste de hiperparámetros para el SVM

# set.seed(1234)
# svm_tune3 <- tune("svm", price_categ1 ~ log_sqft_living + zona + lat, data = datos_train_limpio, kernel = 'radial', 
#                ranges = list(cost = c(0.01, 0.1,1, 10,100), gamma = c(0.01, 0.1, 1, 10, 100)))

load("svm_tune2.RData")

svm_tune <- svm_tune3
svm_tune$best.parameters

##    cost gamma
## 18    1    10

rm(svm_tune3)
best_cost <- svm_tune$best.parameters[1]
best_gamma <- svm_tune$best.parameters[2]

set.seed(1234)
modelo_svm <- e1071::svm(formula = price_categ1 ~ log_sqft_living + zona + lat,
                       data = datos_train_limpio,
                       kernel = "radial",
                       probability =TRUE,
                       gamma = best_gamma$gamma,
                       cost = best_cost$cost)
summary(modelo_svm)

## 
## Call:
## svm(formula = price_categ1 ~ log_sqft_living + zona + lat, data = datos_train_limpio, 
##     kernel = "radial", probability = TRUE, gamma = best_gamma$gamma, 
##     cost = best_cost$cost)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  5340
## 
##  ( 2825 2515 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  B1 C1

6.6.1 Train

tabla_svm_train=table(obs = datos_train_limpio$price_categ1, pred = predict(modelo_svm,datos_train_limpio[,-15], type ="class"))
tabla_svm_train

##     pred
## obs    B1   C1
##   B1 7874  860
##   C1 1002 5383

pred_caras_svm_train=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[1,2])
rec_caras_svm_train=tabla_svm_train[2,2]/(tabla_svm_train[2,2]+tabla_svm_train[2,1])
F_caras_svm_train=(2*pred_caras_svm_train*rec_caras_svm_train)/(pred_caras_svm_train+rec_caras_svm_train)
cat(c('F1 caras: ', F_caras_svm_train), '\n')

## F1 caras:  0.852549889135255

pred_baratas_svm_train=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[2,1])
rec_baratas_svm_train=tabla_svm_train[1,1]/(tabla_svm_train[1,1]+tabla_svm_train[1,2])
F_baratas_svm_train=(2*pred_baratas_svm_train*rec_baratas_svm_train)/(pred_baratas_svm_train+rec_baratas_svm_train)
cat(c('F1 baratas: ', F_baratas_svm_train), '\n')

## F1 baratas:  0.894264622373651

F_svm_train= (F_caras_svm_train+F_baratas_svm_train)/2
cat(c('F1 global: ', F_svm_train), '\n')

## F1 global:  0.873407255754453

accuracy_svm_train = (tabla_svm_train[1,1]+tabla_svm_train[2,2]) / (tabla_svm_train[1,1]+tabla_svm_train[1,2]+tabla_svm_train[2,1]+tabla_svm_train[2,2])
cat(c('Accuracy: ', accuracy_svm_train), '\n')

## Accuracy:  0.876843706594352

En base al ajuste de hiperparámetros, el mejor valor del coste es de 1, y para gamma es de 10. Por último, se evalúa el resultado del modelo, obteniendo una F1-medida y accuracy de 0.87.

7 Evaluación y comparación de modelos

Una vez que se han entrenado y optimizado distintos modelos, se tiene que identificar cuál de ellos consigue mejores resultados para el problema en cuestión, en este caso, predecir si una casa es barata o cara. Con los datos disponibles, existen dos formas de comparar los modelos. Si bien las dos no tienen por qué dar los mismos resultados, son complementarias a la hora de tomar una decisión final.

models_cross = data.frame(
"Modelo"= c('GLM','GAM','KNN','Decision_Tree','Random_Forest','SVM'),
"F_Medida_train" = c(F_reglog_train,F_gam_train,F_knn_train,F_dt_train,F_rf_train,F_svm_train),
"F_Medida_test" = c(F_reglog_test,F_gam_test, F_knn_test,F_dt_test,F_rf_test,F_svm_test),
"F_Medida_validacion" = c(F_reglog_validacion,F_gam_validacion,F_knn_validacion,F_dt_validacion,F_rf_validacion,F_svm_validacion))
models_cross = models_cross[order(-models_cross$F_Medida_train),]
models_cross

##          Modelo F_Medida_train F_Medida_test F_Medida_validacion
## 5 Random_Forest      0.8919936     0.8909831           0.8891082
## 3           KNN      0.8824311     0.8754228           0.8723608
## 6           SVM      0.8734073     0.8620418           0.8739711
## 2           GAM      0.8528492     0.8483874           0.8533047
## 1           GLM      0.8442437     0.8310812           0.8299958
## 4 Decision_Tree      0.8095628     0.8088307           0.8060881

ggplot(data=models_cross, aes(x=reorder(Modelo,F_Medida_train), y=F_Medida_train, fill=Modelo)) +
  xlab("Modelos") +
  geom_bar(stat="identity")

El mejor modelo según la métrica de evaluación considerada (F1-medida para las casas caras y las baratas) es el Random Forest y el peor el Decision Tree, si bien es cierto que todos ellos obtienen un resultado aceptable, entre 0.80 y 0.9.

8 Curva ROC

set.seed(1234)

# REGRESIÓN LOGÍSTICA
predictions_glm <- predict(model_glm, newdata = datos_train_rl, type = "response")
pred_glm <- prediction(predictions_glm, datos_train_rl$price_categ1)
perf_glm <- performance(pred_glm,"spec","sens")

# GAM
predictions_gam <- predict( model_gam_categ, type = "response")
predictions_gam <- as.data.frame(predictions_gam)
pred_gam <- prediction(predictions_gam, datos_train_limpio$price_categ1)
perf_gam <- performance(pred_gam,"spec","sens")

# KNN
predictions_Knn <- knn(scale(datos_train_numeric), scale(datos_train_numeric), cl=datos_train_limpio$price_categ1, k=14, prob = TRUE)
pred_knn <- prediction(attr(predictions_Knn,"prob"), datos_train_limpio$price_categ1)
perf_knn <- performance(pred_knn,"spec","sens")

# ARBOL DE DECISION
predictions_tree <- predict(decisiontree_model, newdata = datos_decision_tree, type = "prob")
pred_tree = prediction(predictions_tree[,2], datos_decision_tree$price_categ1)
perf_tree = performance(pred_tree,"spec","sens")

# RANDOM FOREST
predictions_rf <- predict(randomforest_model, new_data=datos_train_limpio, type = "prob")
pred_rf <- prediction(predictions_rf[,2],datos_train_limpio$price_categ1)
perf_rf <- performance(pred_rf,"spec","sens")

# SVM
predictions_svm <- predict(modelo_svm, newdata=datos_train_limpio, probability = TRUE)
prob_svm<-attr(predictions_svm,"probabilities")
pred_svm <- prediction(prob_svm[,2], datos_train_limpio$price_categ)
perf_svm <- performance(pred_svm,"spec","sens")


plot(perf_glm,col="blue", xlim=c(1,0))
plot(perf_gam,col="black", add = TRUE)
plot(perf_tree, col="red",add = TRUE)
plot(perf_rf,col="green", add = TRUE)
plot(perf_svm, col="yellow",add = TRUE)
plot(perf_knn, col="orange",add = TRUE)
legend(x="right", legend=c("GLM","GAM","DT","RF","SVM","KNN"), fill=c("blue","black","red","green","yellow","orange"), cex=0.8)

9 Conclusiones y líneas futuras

9.1 Conclusiones

Tras el análisis exploratorio y la categorización de la variable respuesta, se realiza un análisis de componentes principales, y se implementan modelos de aprendizaje supervisado y no supervisado. Finalmente, se comparan y evalúan dichos modelos. Adicionalmente, del análisis exploratorio se concluye, que es conveniente categorizar la variable zipcode y reagrupar la variable bathroom.

De los modelos no supervisados (modelos jerárquicos, K-Means y K-medoids), se observa que no hay un patrón definido que permita relacionar de manera clara las características de las casas y la categoría a la que pertenecen según su precio. Con estos métodos no se han obtenido buenos resultados por los motivos que se han comentado.

En cuanto a la métrica que se ha decidido emplear para la evaluación de los métodos supervisados ha sido la F1-medida, por un lado se ha calculado para las casas caras (categoría C1), y por otro lado, para las baratas (B1). Después se ha calculado la media entre ambas. Adicionalmente, se ha hallado la accuracy, que es una métrica de precisión, que tiene en cuenta los aciertos para ambas categorías.

Despues de implementar los modelos supervisados, se concluye que el mejor modelo según la métrica considerada es el Random Forest, con aproximadamente una F1-medida de 0.9, seguido del KNN y el SVM con una F-medida muy similar aproximada de 0.88. En último lugar, se sitúa el árbol de decisión con 0.80.

9.2 Líneas futuras

Implentación de modelos no supervisados basados en similaridad con otras distancias diferentes a la euclidea.
Ajuste de hiperparámetros para algunos modelos para los que no se ha considerado hasta el momento.
Nuevas categorizaciones de price, agrupando por zipcode, considerando a su vez, la agrupación de zipcodes en función de lat y long.
Tener en cuenta algunas variables no consideradas hasta el momento como la fecha de venta de las casas (date).

Machine learning I

Natalia Alonso, Beatriz Visitación y Susana Albarrán

28/04/2020